function [ Ax ] = getAx( i,j,U,V,NJue,NIue,NJut ,h)
%GETAX Returns the value of Ax at the point (i+1/2,j) on the pressure grid.
% [ Ax ] = getAx( i,j,U,V,NJue,NIue,NJut ,h)
%
%   By GDA & CAB
%   Reviewed by GDMT.

% Right and left volumes
U_cen   = U ( Umap (i  ,j  ) );  % center of finite volume
U_right = U ( Umap (i+1,j  ) );
U_left  = U ( Umap (i-1,j  ) );

% Top and bot volumes, considering ghost nodes
if ( j==1 ) U_bot = - U_cen ; else U_bot   = U ( Umap (i  ,j-1) ); end
if (j<=NJue)
    if ( j==NJue && i<=NIue ) U_top   = - U_cen; else U_top   = U ( Umap (i  ,j+1) ); end
else
    if ( j==NJut ) U_top   = - U_cen; else U_top   = U ( Umap (i  ,j+1) ); end
end
% Velocity in Y direction on vertices of the finite volume
V_tr    = V ( Vmap (i  ,j+1) );  % top-right
V_tl    = V ( Vmap (i-1,j+1) );  % top-left
V_br    = V ( Vmap (i  ,j  ) );  % bot-right
V_bl    = V ( Vmap (i-1,j  ) );  % bot-left

Ax  = + ( U_cen + U_right )^2 ...               % Right boundary
    + ( V_tr + V_tl ) * ( U_cen + U_top ) ...   % Top boundary
    - ( U_cen + U_left  )^2 ...                 % Left boundary
    - ( V_br + V_bl ) * ( U_cen + U_bot );      % Bottom boundary

Ax  = Ax/(4*h);


end


